namespace MSL {


// Level 1 Blas Routines


template< typename T >
bool scale(const T & alpha, DistributedDenseMatrix< T > & x) {

      T * message = new T (alpha);
      x.broadcast(message, Muesli::MSL_myId, 1);
      for (std::map<int, Eigen::Matrix< T, Dynamic, Dynamic >*>::const_iterator it = x.submatrices.begin(); 
            it = x.submatrices.end(); it++) {
                  *(it->second) *= *message;
      }
      return true;
      
}


template< typename T >
bool scale(const T & alpha, DistributedDenseVector< T > & x) {

      T * message = new T (alpha);
      x.broadcast(message, Muesli::MSL_myId, 1);
      for (std::map<int, Eigen::Matrix< T, Dynamic, 1 >*>::const_iterator it = x.submatrices.begin(); 
            it = x.submatrices.end(); it++) {
                  *(it->second) *= *message;
      }
      return true;
      
}


template< typename T >
bool scale(const T & alpha, DistributedSparseMatrix< T > & x) {

      T * message = new T (alpha);
      x.broadcast(message, Muesli::MSL_myId, 1);
      for (std::map<int, Eigen::SparseMatrix< T >*>::const_iterator it = x.submatrices.begin(); 
            it = x.submatrices.end(); it++) {
                  *(it->second) *= *message;
      }
      return true;
      
}


template< typename T >
bool scale(const T & alpha, DistributedSparseVector< T > & x) {

      T * message = new T (alpha);
      x.broadcast(message, Muesli::MSL_myId, 1);
      for (std::map<int, Eigen::SparseVector< T >*>::const_iterator it = x.submatrices.begin(); 
            it = x.submatrices.end(); it++) {
                  *(it->second) *= *message;
      }
      return true;
      
}



template< typename T >
bool gaxpy(const T & alpha, const DistributedDenseVector< T > & x, DistributedDenseVector< T > & y) {

      if (x.rowVec) {
            if (y.rowVec) {
                  if (x.m != y.m) {
                        return false;
                  }
            } else {
                  if (x.m != y.n) {
                        return false;
                  }
            }
      } else {
            if (y.rowVec) {
                  if (x.n != y.m) {
                        return false;
                  }
            } else {
                  if (x.n != y.n) {
                        return false;
                  }
            }
      }
            
      for (std::map<int, Eigen::Matrix< T, Dynamic, 1 >*>::const_iterator it = x.submatrices.begin(); 
            it = x.submatrices.end(); it++) {
                   y.submatrices[it->first] += alpha * (*(it->second));
      }
      
      return true;
      
}


template< typename T >
bool gaxpy(const T & alpha, const DistributedDenseVector< T > & x, DistributedSparseVector< T > & y) {

      if (x.rowVec) {
            if (y.rowVec) {
                  if (x.m != y.m) {
                        return false;
                  }
            } else {
                  if (x.m != y.n) {
                        return false;
                  }
            }
      } else {
            if (y.rowVec) {
                  if (x.n != y.m) {
                        return false;
                  }
            } else {
                  if (x.n != y.n) {
                        return false;
                  }
            }
      }
            
      for (std::map<int, Eigen::Matrix< T, Dynamic, 1 >*>::const_iterator it = x.submatrices.begin(); 
            it = x.submatrices.end(); it++) {
                   y.submatrices[it->first] += alpha * (*(it->second));
      }
      
      return true;
      
}


template< typename T >
bool gaxpy(const T & alpha, const DistributedSparseVector< T > & x, DistributedDenseVector< T > & y) {

      if (x.rowVec) {
            if (y.rowVec) {
                  if (x.m != y.m) {
                        return false;
                  }
            } else {
                  if (x.m != y.n) {
                        return false;
                  }
            }
      } else {
            if (y.rowVec) {
                  if (x.n != y.m) {
                        return false;
                  }
            } else {
                  if (x.n != y.n) {
                        return false;
                  }
            }
      }
            
      for (std::map<int, Eigen::SparseMatrix< T >*>::const_iterator it = x.submatrices.begin(); 
            it = x.submatrices.end(); it++) {
                   y.submatrices[it->first] += alpha * (*(it->second));
      }
      
      return true;
      
}


template< typename T >
bool gaxpy(const T & alpha, const DistributedSparseVector< T > & x, DistributedSparseVector< T > & y) {

      if (x.rowVec) {
            if (y.rowVec) {
                  if (x.m != y.m) {
                        return false;
                  }
            } else {
                  if (x.m != y.n) {
                        return false;
                  }
            }
      } else {
            if (y.rowVec) {
                  if (x.n != y.m) {
                        return false;
                  }
            } else {
                  if (x.n != y.n) {
                        return false;
                  }
            }
      }
            
      for (std::map<int, Eigen::SparseMatrix< T >*>::const_iterator it = x.submatrices.begin(); 
            it = x.submatrices.end(); it++) {
                   y.submatrices[it->first] += alpha * (*(it->second));
      }
      
      return true;
      
}


template< T >
T dot(const DistributedDenseVector< T > & x, const DistributedDenseVector< T > & y) {

      if (!(x.distribution->equals(y.distribution)) {
            return x.zero;
      }
      if (x.getElementCount() != other.getElementCount()) {
            return x.zero;
      }

      T localResult = (T) 0.0;
      for (std::map<int, Eigen::Matrix< T, Dynamic, 1 >*>::const_iterator it = x.submatrices.begin(); 
            it = x.submatrices.end(); it++) {
                  localResult += *(it->second).dot(*(y.submatrices[it->first]));
      }      

      localResult = localResult * localResult;
      T globalResult = 0.0;      
      this->allreduce(&localResult, &globalResult, x.add);
      globalResult = std::sqrt(globalResult);
      
      return globalResult;

}


template< T >
T dot(const DistributedDenseVector< T > & x, const DistributedSparseVector< T > & y) {

      if (!(x.distribution->equals(y.distribution)) {
            return x.zero;
      }
      if (x.getElementCount() != other.getElementCount()) {
            return x.zero;
      }

      T localResult = (T) 0.0;
      for (std::map<int, Eigen::Matrix< T, Dynamic, 1 >*>::const_iterator it = x.submatrices.begin(); 
            it = x.submatrices.end(); it++) {
                  localResult += *(it->second).dot(*(y.submatrices[it->first]));
      }      

      localResult = localResult * localResult;
      T globalResult = 0.0;      
      this->allreduce(&localResult, &globalResult, x.add);
      globalResult = std::sqrt(globalResult);
      
      return globalResult;

}


template< T >
T dot(const DistributedSparseVector< T > & x, const DistributedDenseVector< T > & y) {

      if (!(x.distribution->equals(y.distribution)) {
            return x.zero;
      }
      if (x.getElementCount() != other.getElementCount()) {
            return x.zero;
      }

      T localResult = (T) 0.0;
      for (std::map<int, Eigen::SparseMatrix< T >*>::const_iterator it = x.submatrices.begin(); 
            it = x.submatrices.end(); it++) {
                  localResult += *(it->second).dot(*(y.submatrices[it->first]));
      }      

      localResult = localResult * localResult;
      T globalResult = 0.0;      
      this->allreduce(&localResult, &globalResult, x.add);
      globalResult = std::sqrt(globalResult);
      
      return globalResult;

}


template< T >
T dot(const DistributedSparseVector< T > & x, const DistributedSparseVector< T > & y) {

      if (!(x.distribution->equals(y.distribution)) {
            return x.zero;
      }
      if (x.getElementCount() != other.getElementCount()) {
            return x.zero;
      }

      T localResult = (T) 0.0;
      for (std::map<int, Eigen::SparseMatrix< T >*>::const_iterator it = x.submatrices.begin(); 
            it = x.submatrices.end(); it++) {
                  localResult += *(it->second).dot(*(y.submatrices[it->first]));
      }      

      localResult = localResult * localResult;
      T globalResult = 0.0;      
      this->allreduce(&localResult, &globalResult, x.add);
      globalResult = std::sqrt(globalResult);
      
      return globalResult;

}


template< typename T >
T nrm2(const DistributedDenseVector< T > & x) {

      return dot(x, x);
      
}


template< typename T >
T nrm2(const DistributedSparseVector< T > & x) {

      return dot(x, x);
      
}


// Level 2 Blas


// Simple ring-topology MatVec, appropriate for sparse or dense matrices and either dense or sparse vectors;
// "predistributed" parameter specifies whether matrix and vectors have already been appropriately partitioned
// If matrix is not predistributed, MatVec will resdistribute it as a side-effect; in case of sparse matrix, vector 2
// will also be redistributed

template< typename T >
MatVec(const T & alpha,
const DistributedDenseMatrix< T > & M,
const DistributedDenseVector< T > & v1, 
const T & beta,
DistributedDenseVector< T > & v2,
bool predistributed) {

      if (M.getColCount() != v1.getRowCount() || M.getRowCount() != v2.getRowCount()) {
            return;
      }
      
      // Initialize contiguous block-row distribution for M; can assume v1 and v2 already have such a distribution
      if (!(predistributed)) { 
            int np = Muesli::MSL_numOfLocalProcs;
            int n = M.getRowCount();
            int m = M.getColumnCount();
            int r = n / np;
            int c = m;
            int max = n / r;
            if (n % r != 0) {
                  ++max;
            }
            BlockDistribution dist();
            dist.initialize(np, n, m, r, c, max);
            M.blockRedistribute(dist, 0);
      }
            
      // Get shifting distribution for v1
      Distribution * shiftdist;
      
      // Scale v2
      scale(beta, v2);

      // Get pointers to working submatrices; only v1sub changes
      Eigen::Matrix< T, Dymamic, 1 > ** v1sub;
      Eigen::Matrix< T, Dynamic, 1 > ** v2sub;
      v2.getSubmatrices(v2sub);
      Eigen::Matrix< T, Dynamic, Dynamic > ** Msub;
      M.getSubmatrices(Msub);
      
      // Current column ranges
      int startcol, numcols;
      
      // Perform multiplication
      for (int shift = 0; shift < np; ++shift) {
            v1.getSubmatrices(v1sub);
            startcol = v1sub[0]->getRowIndexStart();
            numcols = v1sub[0]->rows();
            *(v2sub[0]) += alpha * Msub[0]->middleCols(startcol, numcols) * *(v1sub[0]);
            shiftdist = v1.oldDist();
            v1.blockRedistribute(*shiftdist, 1);
      }
      
      // Clean up
      delete shiftdist;

}


template< typename T >
MatVec(const T & alpha,
const DistributedDenseMatrix< T > & M,
const DistributedSparseVector< T > & v1, 
const T & beta,
DistributedDenseVector< T > & v2,
bool predistributed) {

      if (M.getColCount() != v1.getRowCount() || M.getRowCount() != v2.getRowCount()) {
            return;
      }
      
      // Initialize contiguous block-row distribution for M; can assume v1 and v2 already have such a distribution
      if (!(predistributed)) { 
            int np = Muesli::MSL_numOfLocalProcs;
            int n = M.getRowCount();
            int m = M.getColumnCount();
            int r = n / np;
            int c = m;
            int max = n / r;
            if (n % r != 0) {
                  ++max;
            }
            BlockDistribution dist();
            dist.initialize(np, n, m, r, c, max);
            M.blockRedistribute(dist, 0);
      }
            
      // Get shifting distribution for v1
      Distribution * shiftdist;
      
      // Scale v2
      scale(beta, v2);

      // Get pointers to working submatrices; only v1sub changes
      Eigen::Matrix< T, Dymamic, 1 > ** v1sub;
      Eigen::Matrix< T, Dynamic, 1 > ** v2sub;
      v2.getSubmatrices(v2sub);
      Eigen::Matrix< T, Dynamic, Dynamic > ** Msub;
      M.getSubmatrices(Msub);
      
      // Current column ranges
      int startcol, numcols;
      
      // Perform multiplication
      for (int shift = 0; shift < np; ++shift) {
            v1.getSubmatrices(v1sub);
            startcol = v1sub[0]->getRowIndexStart();
            numcols = v1sub[0]->rows();
            *(v2sub[0]) += alpha * Msub[0]->middleCols(startcol, numcols) * *(v1sub[0]);
            shiftdist = v1.oldDist();
            v1.blockRedistribute(*shiftdist, 1);
      }
      
      // Clean up
      delete shiftdist;

}


template< typename T >
MatVec(const T & alpha,
const DistributedDenseMatrix< T > & M,
const DistributedDenseVector< T > & v1, 
const T & beta,
DistributedSparseVector< T > & v2,
bool predistributed) {

      if (M.getColCount() != v1.getRowCount() || M.getRowCount() != v2.getRowCount()) {
            return;
      }
      
      // Initialize contiguous block-row distribution for M; can assume v1 and v2 already have such a distribution
      if (!(predistributed)) { 
            int np = Muesli::MSL_numOfLocalProcs;
            int n = M.getRowCount();
            int m = M.getColumnCount();
            int r = n / np;
            int c = m;
            int max = n / r;
            if (n % r != 0) {
                  ++max;
            }
            BlockDistribution dist();
            dist.initialize(np, n, m, r, c, max);
            M.blockRedistribute(dist, 0);
      }    
        
      // Get shifting distribution for v1
      Distribution * shiftdist;
      
      // Scale v2
      scale(beta, v2);

      // Get pointers to working submatrices; only v1sub changes
      Eigen::Matrix< T, Dymamic, 1 > ** v1sub;
      Eigen::Matrix< T, Dynamic, 1 > ** v2sub;
      v2.getSubmatrices(v2sub);
      Eigen::Matrix< T, Dynamic, Dynamic > ** Msub;
      M.getSubmatrices(Msub);
      
      // Current column ranges
      int startcol, numcols;
      
      // Perform multiplication
      for (int shift = 0; shift < np; ++shift) {
            v1.getSubmatrices(v1sub);
            startcol = v1sub[0]->getRowIndexStart();
            numcols = v1sub[0]->rows();
            *(v2sub[0]) += alpha * Msub[0]->middleCols(startcol, numcols) * *(v1sub[0]);
            shiftdist = v1.oldDist();
            v1.blockRedistribute(*shiftdist, 1);
      }
      
      // Clean up
      delete shiftdist;

}


template< typename T >
MatVec(const T & alpha,
const DistributedDenseMatrix< T > & M,
const DistributedSparseVector< T > & v1, 
const T & beta,
DistributedSparseVector< T > & v2,
bool predistributed) {

      if (M.getColCount() != v1.getRowCount() || M.getRowCount() != v2.getRowCount()) {
            return;
      }
      
      // Initialize contiguous block-row distribution for M; can assume v1 and v2 already have such a distribution
      if (!(predistributed)) { 
            int np = Muesli::MSL_numOfLocalProcs;
            int n = M.getRowCount();
            int m = M.getColumnCount();
            int r = n / np;
            int c = m;
            int max = n / r;
            if (n % r != 0) {
                  ++max;
            }
            BlockDistribution dist();
            dist.initialize(np, n, m, r, c, max);
            M.blockRedistribute(dist, 0);
      }
      
      // Get shifting distribution for v1
      Distribution * shiftdist;
      
      // Scale v2
      scale(beta, v2);

      // Get pointers to working submatrices; only v1sub changes
      Eigen::Matrix< T, Dymamic, 1 > ** v1sub;
      Eigen::Matrix< T, Dynamic, 1 > ** v2sub;
      v2.getSubmatrices(v2sub);
      Eigen::Matrix< T, Dynamic, Dynamic > ** Msub;
      M.getSubmatrices(Msub);
      
      // Current column ranges
      int startcol, numcols;
      
      // Perform multiplication
      for (int shift = 0; shift < np; ++shift) {
            v1.getSubmatrices(v1sub);
            startcol = v1sub[0]->getRowIndexStart();
            numcols = v1sub[0]->rows();
            *(v2sub[0]) += alpha * Msub[0]->middleCols(startcol, numcols) * *(v1sub[0]);
            shiftdist = v1.oldDist();
            v1.blockRedistribute(*shiftdist, 1);
      }
      
      // Clean up
      delete shiftdist;

}


// Redistribution for sparse matrices is different, uses load balancing
template< typename T >
MatVec(const T & alpha,
const DistributedSparseMatrix< T > & M,
const DistributedDenseVector< T > & v1, 
const T & beta,
DistributedDenseVector< T > & v2,
bool predistributed) {

      if (M.getColCount() != v1.getRowCount() || M.getRowCount() != v2.getRowCount()) {
            return;
      }
      
      // Redistribute for sparse matvec if necessary
      if (!(predistributed)) {
            M.partitionForRowwiseBalance();
            std::vector< int > newcuts = M.shareCuts();
            v2.rowRedistribute(newcuts);
            v1.partitionForRowwiseBalance();
      }
      
      // Get shifting distribution for v1
      Distribution * shiftdist;
      
      // Scale v2
      scale(beta, v2);

      // Get pointers to working submatrices; only v1sub changes
      Eigen::Matrix< T, Dymamic, 1 > ** v1sub;
      Eigen::Matrix< T, Dynamic, 1 > ** v2sub;
      v2.getSubmatrices(v2sub);
      Eigen::SparseMatrix< T > ** Msub;
      M.getSubmatrices(Msub);
      
      // Current column ranges
      int startcol, numcols;
      
      // Perform multiplication
      for (int shift = 0; shift < np; ++shift) {
            v1.getSubmatrices(v1sub);
            startcol = v1sub[0]->getRowIndexStart();
            numcols = v1sub[0]->rows();
            *(v2sub[0]) += alpha * Msub[0]->middleCols(startcol, numcols) * *(v1sub[0]);
            shiftdist = v1.oldDist();
            v1.blockRedistribute(*shiftdist, 1);
      }
      
      // Clean up
      delete shiftdist;

}


template< typename T >
MatVec(const T & alpha,
const DistributedSparseMatrix< T > & M,
const DistributedSparseVector< T > & v1, 
const T & beta,
DistributedDenseVector< T > & v2) {

      if (M.getColCount() != v1.getRowCount() || M.getRowCount() != v2.getRowCount()) {
            return;
      }
      
      // Redistribute for sparse matvec if necessary
      if (!(predistributed)) {
            M.partitionForRowwiseBalance();
            std::vector< int > newcuts = M.shareCuts();
            v2.rowRedistribute(newcuts);
            v1.partitionForRowwiseBalance();
      }
      
      // Get shifting distribution for v1
      Distribution * shiftdist;
      
      // Scale v2
      scale(beta, v2);

      // Get pointers to working submatrices; only v1sub changes
      Eigen::Matrix< T, Dymamic, 1 > ** v1sub;
      Eigen::Matrix< T, Dynamic, 1 > ** v2sub;
      v2.getSubmatrices(v2sub);
      Eigen::SparseMatrix< T > ** Msub;
      M.getSubmatrices(Msub);
      
      // Current column ranges
      int startcol, numcols;
      
      // Perform multiplication
      for (int shift = 0; shift < np; ++shift) {
            v1.getSubmatrices(v1sub);
            startcol = v1sub[0]->getRowIndexStart();
            numcols = v1sub[0]->rows();
            *(v2sub[0]) += alpha * Msub[0]->middleCols(startcol, numcols) * *(v1sub[0]);
            shiftdist = v1.oldDist();
            v1.blockRedistribute(*shiftdist, 1);
      }
      
      // Clean up
      delete shiftdist;

}


template< typename T >
MatVec(const T & alpha,
const DistributedSparseMatrix< T > & M,
const DistributedDenseVector< T > & v1, 
const T & beta,
DistributedSparseVector< T > & v2) {

      if (M.getColCount() != v1.getRowCount() || M.getRowCount() != v2.getRowCount()) {
            return;
      }
      
      // Redistribute for sparse matvec if necessary
      if (!(predistributed)) {
            M.partitionForRowwiseBalance();
            std::vector< int > newcuts = M.shareCuts();
            v2.rowRedistribute(newcuts);
            v1.partitionForRowwiseBalance();
      }
      
      // Get shifting distribution for v1
      Distribution * shiftdist;
      
      // Scale v2
      scale(beta, v2);

      // Get pointers to working submatrices; only v1sub changes
      Eigen::Matrix< T, Dymamic, 1 > ** v1sub;
      Eigen::Matrix< T, Dynamic, 1 > ** v2sub;
      v2.getSubmatrices(v2sub);
      Eigen::SparseMatrix< T > ** Msub;
      M.getSubmatrices(Msub);
      
      // Current column ranges
      int startcol, numcols;
      
      // Perform multiplication
      for (int shift = 0; shift < np; ++shift) {
            v1.getSubmatrices(v1sub);
            startcol = v1sub[0]->getRowIndexStart();
            numcols = v1sub[0]->rows();
            *(v2sub[0]) += alpha * Msub[0]->middleCols(startcol, numcols) * *(v1sub[0]);
            shiftdist = v1.oldDist();
            v1.blockRedistribute(*shiftdist, 1);
      }
      
      // Clean up
      delete shiftdist;

}


template< typename T >
MatVec(const T & alpha,
const DistributedSparseMatrix< T > & M,
const DistributedSparseVector< T > & v1, 
const T & beta,
DistributedSparseVector< T > & v2) {

      if (M.getColCount() != v1.getRowCount() || M.getRowCount() != v2.getRowCount()) {
            return;
      }
      
      // Redistribute for sparse matvec if necessary
      if (!(predistributed)) {
            M.partitionForRowwiseBalance();
            std::vector< int > newcuts = M.shareCuts();
            v2.rowRedistribute(newcuts);
            v1.partitionForRowwiseBalance();
      }
      
      // Get shifting distribution for v1
      Distribution * shiftdist;
      
      // Scale v2
      scale(beta, v2);

      // Get pointers to working submatrices; only v1sub changes
      Eigen::Matrix< T, Dymamic, 1 > ** v1sub;
      Eigen::Matrix< T, Dynamic, 1 > ** v2sub;
      v2.getSubmatrices(v2sub);
      Eigen::SparseMatrix< T > ** Msub;
      M.getSubmatrices(Msub);
      
      // Current column ranges
      int startcol, numcols;
      
      // Perform multiplication
      for (int shift = 0; shift < np; ++shift) {
            v1.getSubmatrices(v1sub);
            startcol = v1sub[0]->getRowIndexStart();
            numcols = v1sub[0]->rows();
            *(v2sub[0]) += alpha * Msub[0]->middleCols(startcol, numcols) * *(v1sub[0]);
            shiftdist = v1.oldDist();
            v1.blockRedistribute(*shiftdist, 1);
      }
      
      // Clean up
      delete shiftdist;

}


